% To do likelihood ratio test
clear all;
clc;
disp('To do likelihood ratio test_MK4');

%load MK2 data (Both tests based on same data)
% load PathCountTotal_MK4_TP06YRJfinR5G1000kMR7.mat;
% load PathCountAvgProb_MK4_TP06YRJfinR5G1000kMR7.mat;
load ('MK4_TP06YRJfinR5G2501000kMR7.mat','PathCountTotal','PathCountAvgProb');
filename='TP06_Y_RJ_fin_R5_G250_1000k_MR7_Likelihood Ratio Test_MK3_MK4.xls';

% take the avg data from PathCountAvgProb (MK4).
N4(:,:,:,:,:) = PathCountTotal(1,:,:,:,:,:);
P4(:,:,:,:,:) = PathCountAvgProb(1,:,:,:,:,:);
ttlData = sum(N4(:));

% to calculate log (LMK3) hyphothesis
N3(:,:,:,:) = sum(N4(1:7,:,:,:,:));
P3=zeros(7,7,7,8);
LLMK3=0;
for i=1:7
    for j=1:7
        for k=1:7
            for d=1:8
                P3(i,j,k,d)=N3(i,j,k,d)./sum(N3(i,j,k,1:8));
                if N3(i,j,k,d)~=0
                    Y1 = N3(i,j,k,d);
                    Y2 = P3(i,j,k,d);
                    Y =log(Y2).* Y1;
%                     Y =log(P3(i,j,k,d)).* N3(i,j,k,d);
                    LLMK3=LLMK3+Y;
                else
                    LLMK3=LLMK3+0;
                end
            end
        end
    end
end


% to calculate log (LMK4) hyphothesis
LLMK4=0;
for i=1:7
    for j=1:7
        for k=1:7
            for d=1:7
                for t=1:8
                    if N4(i,j,k,d,t)~=0
                        YY =log(P4(i,j,k,d,t)).* N4(i,j,k,d,t);
                        LLMK4=LLMK4+YY;
                    else
                        LLMK4=LLMK4+0;
                    end
                end
            end
        end
    end
end


%Likelihood Ratio Test
invalid_Num_N4 = prod(size(find(N4==0)));
invalid_Num_N3 = prod(size(find(N3==0)));
Num_N3=prod(size(P3));
Num_N4=prod(size(P4));
%dof = Num_N4 - Num_N3 - invalid_Num_N3 - invalid_Num_N4;
dof = (Num_N4 - invalid_Num_N4) - (Num_N3 - invalid_Num_N3);
[LRh,LRp,LRstat,cV] = lratiotest(LLMK4,LLMK3,dof);


% PP2(:,:)=P2(2,:,:);
% NN2(:,:)=N2(2,:,:); 
% PP3(:,:)=P3(2,2,:,:);
% NN3(:,:)=N3(2,2,:,:); 
% sheetlist1=strcat('Prob of ijk_pathR2');
% sheetlist2=strcat('Prob of ijkd_pathR2R2');
% sheetlist3=strcat('ObservationNum N2_pathR2');
% sheetlist4=strcat('ObservationNum N3_pathR2R2');
sheetlist5=strcat('Parameters');
% xlswrite(filename, PP2, sheetlist1);
% xlswrite(filename, PP3, sheetlist2);
% xlswrite(filename, NN2, sheetlist3);
% xlswrite(filename, NN3, sheetlist4);
d = {'LRh','LRp','LRstat','cV','dof','LLMK3','LLMK4','ttlData'; LRh, LRp, LRstat,cV,dof,LLMK3,LLMK4,ttlData};
xlswrite(filename, d, sheetlist5)


